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I. INTRODUCTION 

If you were given only one method for electronic struc- 
ture calculations, which method would you choose? It 
certainly depends on your field of research, but, on aver- 
age, the usual choice is probably density-functional the- 
ory (DFT). It is not the best method for everything, 
but its efficiency and accuracy are suitable for most pur- 
poses. You may disagree with this argument, but already 
DFT community's size is a convincing evidence of DFT's 
importance — it is even among the few quantum mechan- 
ical methods used in industry. 

Things were different before. When computational re- 
sources were modest and DFT functionals inaccurate, 
classical force fields, scmicmpirical tight-binding, and jel- 
lium DFT calculations were used. Today tight-binding is 
mostly familiar from solid-state textbooks as a method 
for modeling band-structures, with one to several fitted 
hopping parameters^. However, tight-binding could be 
used better than this more often even today especially 
as a method to calculate total energies. Particularly 
useful for total energy calculations is density-functional 
tight-binding, which is parametrized directly using DFT, 
and is hence rooted in first principles deeper than other 
tight-binding flavors. It just happened that density- 
functional tight-binding came into existence some ten 
years ago 2 ^ 3 - when atomistic DFT calculations for real- 
istic system sizes were already possible. With DFT as a 
competitor, the DFTB community never grew large. 

Despite being superseded by DFT, DFTB is still use- 
ful in many ways: i) In calculations of large systems*^ 
Computational scaling of DFT limits system sizes, while 
better scaling is easier to achieve with DFTB. ii) Ac- 
cessing longer time scales. Systems that are limited for 
optimization in DFT can be used for extensive stud- 
ies on dynamical properties in DFTB 6 . iii) Structure 
search and general trends 7 -. Where DFT is limited to 
only a few systems, DFTB can be used to gather statis- 
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tics and trends from structural families. It can be used 
also for pre-screening of systems for subsequent DFT 
calculations^, iv) Method development. The formalism 
is akin to that of DFT, so methodology improvements, 
quick to test in DFTB, can be easily exported and ex- 
tended into DFT—. v) Testing, playing around, learn- 
ing and teaching. DFTB can be used simply for playing 
around, getting the feeling of the motion of atoms at a 
given temperature, and looking at the chemical bonding 
or realistic molecular wavefunctions, even with real-time 
simulations in a classroom — DFTB runs easily on a lap- 
top computer. 

DFTB is evidently not an ab initio method since it 
contains parameters, even though most of them have 
a theoretically solid basis. With parameters in the 
right place, however, computational effort can be re- 
duced enormously while maintaining a reasonable accu- 
racy. This is why DFTB compares well with full DFT 
with minimal basis, for instance. Scmicmpirical tight- 
binding can be accurately fitted for a given test set, 
but transferability is usually worse; for general purposes 
DFTB is a good choice among the different tight-binding 
flavors. 

Despite having its origin in DFT, one has to bear in 
mind that DFTB is still a tight-binding method, and 
should not generally be considered to have the accu- 
racy of full DFT. Absolute transferability can never be 
achieved, as the fundamental starting point is tightly 
bound electrons, with interactions ultimately treated per- 
turbativcly. DFTB is hence ideally suited for covalent 
systems such as hydrocarbons<2iii Nevertheless, it does 
perform surprisingly well describing also metallic bond- 
ing with delocalized valence electron s] 8 ' 12 . 

This article is neither a historical review of different fla- 
vors of tight-binding, nor a review of even DFTB and its 
successes or failures. For these purposes there are several 
good reviews, see, for example Refs. [Mlllill Some 
ideas were around beforei 7 -^^, but the DFTB method 
in its present formulation, presented also in this article, 
was developed in the mid-QO's 2 - ^ 20 ' 21 ' 22 ' 23 . The main 
architects behind the machinery were Eschrig, Seifert, 
Frauenheim, and their co-workers. We apologize for 
omitting other contributors — consult Ref. |l3| for a more 
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organized literature review on DFTB. 

Instead of reviewing, we intend to present DFTB in a 
pedagogical fashion. By occasionally being more explicit 
than usual, we derive the approximate formalism from 
DFT in a systematic manner, with modest referencing 
to how the same approximations were done before. The 
approach is practical: we present how to turn DFT into 
a working tight-binding scheme where parametrizations 
are obtained from well-defined procedures, to yield actual 
numbers — without omitting ugly details hiding behind 
the scenes. Only basic quantum mechanics and selected 
concepts from density-functional theory are required as 
prc-requisitc. 

The DFTB parametrization process is usually pre- 
sented as superficially easy, while actually it is difficult, 
especially regarding the fitting of short-range repulsion 
(Sec. IIV|) . In this article we want to present a systematic 
scheme for fitting the repulsion, in order to accurately 
document the way the parametrization is done. 

We discuss the details behind tight-binding formalism, 
like Slater-Koster integrals and transformations, largely 
unfamiliar for density-functional community; some read- 
ers may prefer to skip these detailed appendices. Because 
one of DFTB's strengths is the transparent electronic 
structure, in the end we also present selected analysis 
tools. 

We concentrate on ground-state DFTB, leaving time- 
dependen t 24 ' 25 ! 26 ' 27 or linear response 2 ^ formalisms out- 
side the discussion. We do not include spin in the for- 
malism. Our philosophy lies in the limited benefits of im- 
proving upon spin-paired self-consistent -charge DFTB. 
Admittedly, one can adjust parametrizations for certain 
systems, but the tight-binding formalism, especially the 
presence of the repulsive potential, contains so many ap- 
proximations that the next level in accuracy, in our opin- 
ion, is full DFT. 

This philosophy underlies hotbit^ software. It is an 
open-source DFTB package, released under the terms of 
GNU general public license^. It has an interface with the 
atomic simulation environment (ASE)^ 1 -, a python mod- 
ule for multi-purpose atomistic simulations. The ASE 
interface enables simulations with different levels of the- 
ory, including many DFT codes or classical potentials, 
with DFTB being the lowest-level quantum-mechanical 
method, hotbit is built upon the theoretical basis de- 
scribed here — but we avoid technical issues related prac- 
tical implementations having no scientific relevance. 



II. THE ORIGINS OF DFTB 
A. Warm-up 

We begin by commenting on practical matters. The 
equations use h 2 /m e = 47re = e = 1. This gives Bohr 
radius as the unit of length (ag = 0.5292 A) and Hartree 
as the unit of energy (Ha= 27.2114 eV). Selecting atomic 
mass unit (u — 1.6605 • 10~ 27 kg) the unit of mass, the 



unit of time becomes 1.0327 fs, appropriate for molecu- 
lar dynamics simulations. Some useful fundamental con- 
stants are then h = yjm e /u = 0.0234, k B = 3.1668 -lO" 6 , 
and £o = l/(47r), for instance. 

Electronic eigenstates are denoted by ip a , and (pseudo- 
atomic) basis states (p^, occasionally adopting Dirac's no- 
tation. Greek letters [i, v are indices for basis states, 
while capital Roman letters /, J are indices for atoms; 
notation [i £ I stands for orbital fi that belongs to atom 
/. Capital R denotes nuclear positions, with the posi- 
tion of atom / at Rj, and displacements Rjj = \Rij\ = 
\Rj — Rj\. Unit vectors are denoted by R = R/\R\. 

In other parts our notation remains conventional; de- 
viations are mentioned or made self-explanatory. 

B. Starting Point: Full DFT 

The derivation of DFTB from DFT has been presented 
several times; see, for example Refs. [l3lfl8l and H. We 
do not want to be redundant, but for completeness we 
derive the equations briefly; our emphasis is on the final 
expressions. 

We start from the total energy expression of interacting 
electron system 

E = T + E cxt + E cc + E u , (1) 

where T is the kinetic energy, E cxt the external in- 
teraction (including electron-ion interactions), E ee the 
electron-electron interaction, and En ion-ion interaction 
energy. Here En contains terms like ZjZ v j/\Rj — Rj\, 
where Zf is the valence of the atom /, and other con- 
tributions from the core electrons. In density- functional 
theory the energy is a functional of the electron density 
n(r), and for Kohn-Sham system of non- interacting elec- 
trons the energy can be written as 

E[n{r)} =T S + E cxt + E H + E xc + E n , (2) 

where T s is the non-interacting kinetic energy, Eh is the 
Hartree energy, and E xc = (T — T s ) + (E ee — E H ) is the 
exchange-correlation (xc) energy, hiding all the difficult 
many-body effects. More explicitly, 

E[n] =5><*.l (4 V2 + V <* + \ J l^^f) 
+ E xc [n]+E n , 

(3) 

where f a £ [0, 2] is the occupation of a single-particle 
state ip a with energy e a , usually taken from the Fermi- 
function (with factor 2 for spin) 

f a = f(e a ) = 2 ■ [exp(e a - fx)/k B T + l]" 1 (4) 

with chemical potential fi chosen such that f a = num- 
ber of electrons. The Hartree potential 

M»](r) = /p^|, (5) 
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is a classical electrostatic potential from given n(r); for 
brevity we will use the notation J d 3 r — > J , J d 3 r' — > J , 
n(r) — > n, and n(r') — ► n'. With this notation the Kohn- 
Sham DFT energy is, once more, 



em =X)/ <v«i (-^ v2 + / 



) m 



(6) 



r — r' 



So far everything is exact, but now we start approx- 
imating. Consider a system with density rio(r) that is 
composed of atomic densities, as if atoms in the system 
were free and neutral. Hence no(r) contains (artificially) 
no charge transfer. The density no(r) does not minimize 
the functional E[n(r)], but neighbors the true minimiz- 
ing density n min {r) = n (r) + Sn (r), where 6n (r) is 
supposed to be small. Expanding E[n] at no(r) to second 
order in fluctuation Sn(r) the energy reads 

E[Sn] « J2 fM - \ v ' 2 + v ^ + V h[M + V xc [noUa) 



5 2 E xc [n ] 1 



8ndn' 



SnSn' 



\ J V H [no}(r)n {r) + E xc [n ] + E n 



V xc [n ](r)n (r), 



(7) 



while linear terms in Sn vanish. The first line in Eq. ([7]) 
is the band-structure energy 



no \\ipa), 



(8) 



where the Hamiltonian H° = H[no] itself contains no 
charge transfer. The second line in Eq. ([7]) is the energy 
from charge fluctuations, being mainly Coulomb interac- 
tion but containing also xc-contributions 

The third and fourth lines in Eq. ([7J are collectively 
called the repulsive energy 



E lcp = - i J V H [n }{r)n (r) + E xc [n ] + E u 



(10) 



V xc [n }(r)n (r), 



because of the ion-ion repulsion term. Using this termi- 
nology the energy is 



E[Sn] = E BS [Sn] + ^coui[^n] + E lep . 



(11) 



Before switching into tight-binding description, we dis- 
cuss E rcp and -E C oui separately and introduce the main 
approximations. 



C. Repulsive Energy Term 

In Eq. (|10p we lumped four terms together and referred 
them to as repulsive interaction. It contains the ion-ion 
interaction so it is repulsive (at least at small atomic 
distances), but it contains also xc-interactions, so it is 
a complicated object. At this point we adopt manners 
from DFT: we sweep the most difficult physics under the 
carpet. You may consider E Tep as practical equivalent to 
an ccc-functional in DFT because it hides the cumbersome 
physics, while we approximate it with simple functions. 

For example, consider the total volumes in the first 
term, the Hartree term 



(12) 



divided into atomic volumes; the integral becomes a sum 
over atom pairs with terms depending on atomic num- 
bers alone, since no(r) depends on them. We can hence 
approximate it as a sum of terms over atom pairs, where 
each term depends only on elements and their distance, 
because no(r) is spherically symmetric for free atoms. 
Similarly ion-ion repulsions 



Z}Z v j 



\Rj-Ri\ 



Ru 



(13) 



depend only on atomic numbers via their valence num- 
bers Z\ . Using similar reasoning for the remaining terms 
the repulsive energy can be approximated as 



(14) 



KJ 



For each pair of atoms I J we have a repulsive func- 
tion U r gp(i?) depending only on atomic numbers. Note 
that E-ccp contains also on-site contributions, not only the 
atoms' pair-wise interactions, but these depend only on 
no(r) and shift the total energy by a constant. 

The pair- wise repulsive functions V^ ep {R) are obtained 
by fitting to high-level theoretical calculations; detailed 
description of the fitting process is discussed in IIVI 



D. Charge Fluctuation Term 

Let us first make a side road to recall some general 
concepts from atomic physics. Generally, the atom en- 
ergy can be expressed as a function of Aq extra electrons 
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\dAq 



2 \dAq 2 



(15) 



= E - X Aq+-UAq 2 . 



The (negative) slope of E(Aq) at Aq = is given by 
the (positive) electronegativity, which is usually approx- 
imated as 



X w (IE + EA)/2, 



(16) 
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where IE the ionization energy and EA the electron 
affinity. The (upward) curvature of E(Aq) is given by 
the Hubbard U, which is 



Ufa IE- EA, 



(17) 



and is twice the atom absolute hardness rj (U = 2r])21. 
Electronegativity comes mainly from orbital energies rel- 
ative to the vacuum level, while curvature effects come 
mainly from Coulomb interactions. 

Let us now return from our side road. The energy in 
Eq. ([9|) comes from Coulomb and xc-interactions due to 
fluctuations Sn(r), and involves a double integrals over 
all space. Consider the space V divided into volumes Vi 
related to atoms /, such that 



Vi = V and 



E 



Vi 



(18) 



We never precisely define what these volumes V/ exactly 
are — they are always used qualitatively, and the usage 
is case-specific. For example, volumes can be used to 
calculate the extra electron population on atom / as 



A qi 



6n(r)d 3 



(19) 



By using these populations we can decompose 5n into 
atomic contributions 



Sn(r) = 2_] AqiSrii(r), 



i 



(20) 



such that each Snj(r) is normalized, f v Sni(r)d 3 r = 1. 



Note that Eqs. (| 1Q[) and ([20]) are internally consistent. 
Ultimately, this division is used to convert the double 
integral in Eq. ^ into sum over atoms pairs I J, and 
integrations over volumes J v J v ■ 
First, terms with I = J are 



Mi 



S 2 E xc [n ] 
SnSn' 



1 



r — r 



SmSnj. (21) 



Term depends quadratically on Aqj and by comparing 
to Eq. (|15p we can see that integral can be approximated 
by U. Hence, terms with I = J become \U\Aq 2 . 

Second, when I =/= J rrc-contributions will vanish for 
local xc-functionals for which 
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oc 5{r 



Sn(r)Sn(r') 
and the interaction is only electrostatic, 

' SnjSn'j 



1 



AqiAqj 



v f J v.. 



r — r 



(22) 



(23) 



between extra atomic populations Aqj and Aqj . Strictly 
speaking, we do not know what the functions Snj(r) are. 
However, assuming spherical symmetry, they tell how 
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FIG. 1: (color online) (a) The change in density profile for C 
upon charging. Atom is slightly charged (—2/3 e and +1 e), 
and we plot the averaged Sn(r) = |n ± (r) —no(r)\ (shadowed), 
where n ± (r) is the radial electron density for charged atom, 
and no(r) is the radial electron density for neutral atom. This 
is compared to the Gaussian profile of Eq. (|24p with FWHM= 
1.329/C/ where U is given by Eq. (|17|) , The change in density 
near the core is irregular, but the behavior is smooth for up to 
~ 90 % of the density change, (b) The interaction energy of 
two spherically symmetric Gaussian charge distributions with 
equal FWHM, =FWHMj = 1.329/17 with U = 1, as given by 
Eq. (J26J). With R u > FWHM/ interaction is Coulomb-like, 
and approaches U as Ru — ► 0. 



the density profile of a given atom changes upon charg- 
ing. By assuming functional form for the profiles Sni(r), 
the integrals can be evaluated. We choose a Gaussian 
profile^, 



5ni(r) 



1 



(271^)3/: 



■ exp 



2a] 



where 



07 = 



FWHM/ 



V8hT2 



(24) 



(25) 



and FWHM/ is the full width at half maximum for the 
profile. This choice of profile is justified for a carbon atom 
in Fig. QJi. With these assumptions, the Coulomb energy 
of two spherically symmetric Gaussian charge distribu- 
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tions in Eq. (|23[) can be calculated analytically to yield 



-=jij(Rij), (26) 



v Jv 



where 



r — r 



Cu = 



41n2 



FWHMj + FWHM 2 j ' 



(27) 



In definition Q26|) we integrate over whole spaces because 
dni's are strongly localized. The function (j2"6")) is plotted 
in Fig. [TJd, where we see thata when R ^> FWHM we 
get point-like l/i?-interaction. Furthermore, as R — > 0, 
7 — > C • 2/y / 7r, which gives us a connection to on-site 
interactions: if J = J 



777 (i? 7/ = 0) = 



!ln2 



1 



7T FWHM/ 



(28) 



This is the on-site Coulomb energy of extra population 
on atom I. Comparing to the 1 = 3 case above, we can 
approximate 



FWHM, = 



81n2 1 



1.329 



Vj Ui 



(29) 



We interpret Eq. (|28[) as: narrower atomic charge distri- 
butions causes larger costs to add or remove electrons; 
for a point charge charging energy diverges, as it should. 

Hence, from absolute hardness, by assuming only Cou- 
lombic origin, we can estimate the sizes of the charge 
distributions, and these sizes can be used to estimate 
Coulomb interactions also between the atoms. U and 
FWHM are coupled by Eq. ([2"5jl. and hence for each el- 
ement a single parameter Ui, which can be found from 
standard tables, determines all charge transfer energet- 
ics. 

To conclude this subsection, the charge fluctuation in- 
teractions can be written as 



Seoul = ^ ^ 7ij(Rij)AqiAqj, 



i j 



where 



Jij(Rij) 




(30) 



(31) 



E. TB Formalism 

So far the discussion has been without references to 
tight-binding description. Things like eigenstates \ip a ) or 
populations Aqi can be understood, but what are they 
exactly? 

As mentioned above, we consider only valence elec- 
trons; the repulsive energy contains all the core electron 



effects. Since tight-binding assumes tightly bound elec- 
trons, we use minimal local basis to expand 



(32) 



Minimality means having only one radial function for 
each angular momentum state: one for s-states, three for 
p-states, five for d-states, and so on. We use real spheri- 
cal harmonics familiar from chemistry — for completeness 
they are listed in Table ITTl in Appendix lAl 

With this expansion the band-structure energy be- 
comes 



where 



e bs = E $ a E c m* c °^°^ 



(33) 



(34) 



The tight-binding formalism is adopted by accepting the 
matrix elements H® themselves as the principal param- 
eters of the method. This means that in tight-binding 
spirit the matrix elements H® u are just numbers. Calcu- 
lation of these matrix elements is discussed in Section lHll 
with details left in Appendix [Cl 

How about the atomic populations Aqi? Using the 
localized basis, the total number of electrons on atom I 



is 
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a •' V i 



(35) 



If neither /i nor v belong to /, the integral is roughly 
zero, and if both \x and v belong to /, the integral is 
approximately since orbitals on the same atom are 
orthonormal. If [i belongs to / and v to some other atom 
J, the integral becomes 



1 



Vi 



1 



WtevW ~ 2 / ^W^W = 2^' (36) 



as suggested by Fig. [2l where £ M „ = (<Pfi\<p v ) is the over- 
lap of orbitals [i and v. Charge on atom I is 



(37) 



where c.c. stands for complex conjugate. Hence Aqi = 
qi — qj, where q° is the number of valence electrons for 
a neutral atom. This approach is called the Mulliken 
population analysis^. 

Now were are ready for the final energy expression, 



\ E lu(Rij)A qi Aqj + E V X %{R U ) 



(38) 



u 



KJ 



6 




Vi 

FIG. 2: Integrating the overlap of local orbitals. The large 
shaded area denotes the volume Vi of atom /, spheres rep- 
resent schematically the spatial extent of orbitals, and the 
small hatched area denotes the overlap region. With v 6 J 
and [i £ /, the integration of <^ M (r)* tp v (r) over atom Fs vol- 
ume Vi misses half of the overlap, since the other half is left 
approximately to Vj. 

where everything is, in principle, defined. We find 
the minimum of this expression by variation of 8{E — 
Tlia £ a(tpa\' l Pa)), where e a arc undetermined Lagrange 
multipliers, constraining the wave function norms, and 
obtain 

J2^(H^-e a S flu ) = 0, (39) 

for all a and \x. This equation for the coefficients is 
the Kolm-Sham equation -equivalent in DFTB. Here 

Hf+u = H° v + -S^v '^(■Jik + 7jk)^Qk, i" S I v e J. 

K 

(40) 

By noting that the electrostatic potential on atom / due 
to charge fluctuations is e/ = ^2 K r yiK^-QK, the equation 
above can be written as 

= H® v + h^S^v, (41) 

where 

hl„ = \{ei + ej), fi£l veJ. (42) 

This expression suggests a reasonable interpretation: 
charge fluctuations shift the matrix element accord- 
ing to the averaged electrostatic potentials around or- 
bitals (i and v. As in Kohn-Sham equations in DFT, also 
Eqs. (|39|) and ([40|) have to be solved self-consistently: 
from a given initial guess for {Ag/} one obtains /i* and 
H^, then by solving Eq. (f39]) one obtains new {c°}, 
and, finally, new {Ag/}, iterating until self-consistency 
is achieved. The number of iterations required for con- 
vergence is usually markedly less than in DFT, albeit 
similar convergence problems are shared. 

Atomic forces can be obtained directly by taking gradi- 
ents of Eq. (J5HJ) with respect to coordinates (parameters) 
Ri. We get (with Vj = d/dRj) 

-A qi ^(V /7/ j)Ag J - V/£: rep , (43) 



where the gradients of "fu arc obtained analytically from 
Eq. (|2l)|) , and the gradients of and are obtained 
numerically from an interpolation, as discussed in Ap- 
pendix [cj 

III. MATRIX ELEMENTS 

Now we discuss how to calculate the matrix elements 
H^ v and S^. In the main text we describe only main 
ideas; more detailed issues are left to Appendices El IB1 
and [C] 

A. The Pseudo-atom 

The minimal basis functions ip^ in the expansion (1321) . 

W(r') - R^{r)Y^e^) (r' = Rj + r, // G I), (44) 

with real spherical functions Y^(6 1 (p) as defined in Ap- 
pendix should robustly represent bound electrons in 
a solid or molecule, which is what we ultimately want to 
simulate. Therefore orbitals should not come from free 
atoms, as they would be too diffuse. To this end, we 
use the orbitals from a pseudo-atom, where an additional 
confinement potential V r con f(?*) is added to the Hamilto- 
nian 

_ I V 2 - - + V H {r) + V xc (r) + V coni (r). (45) 
2 r 

This additional, spherically symmetric confinement cuts 
the orbitals' diffuse tails off and makes a compact basis — 
and ultimately better basis^i — for the wave function ex- 
pansion. 

A general, spherically symmetric environment can be 
represented by a potential 

oo 

V con{ (r) = ^v 2i r 2i , (46) 

i=0 

where the odd terms disappear because the potential has 
to be smooth at r = 0. Since the first vq term is just a 
constant shift, the first non-trivial term is v^r 2 . To first 
approximation we hence choose the confining potential 
to be of the form 

V co ^{r)=(£\ , (47) 

where is a parameter. The quadratic form for the 
confinement has appeared before^, but also other forms 
have been analyzed^. While different forms can be con- 
sidered for practical reasons, they have only little effect 
on DFTB performance. The adjustment of the parame- 
ter r is discussed in Section IVl 

The pseudo-atom is calculated with DFT only once for 
a given confining potential. This way we get (p^ 's (more 
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precisely, i? M 's), the localized basis functions, for later 
use in matrix element calculations. 

One technical detail we want to point out here con- 
cerns orbital conventions. Namely, once the orbitals tp^ 
are calculated, their sign and other conventions should 
never change. Fixed convention should be used in all 
simultaneously used Slater-Koster tables; using different 
conventions for same elements gives inconsistent tables 
that are plain nonsense. The details of our conventions, 
along with other technical details of the pseudo-atom cal- 
culations, are discussed in Appendix [A"l 



B. Overlap Matrix Elements 

Using the orbitals from pseudo-atom calculations, we 
need to calculate the overlap matrix elements 

-V = y^(r)V,(r)d 3 r. (48) 

Since orbitals are chosen real, the overlap matrix is real 
and symmetric. 

The integral with ip^ at Rj and <p v at Rj can be calcu- 
lated also with ip^ at the origin and tp v at Rij. Overlap 
will hence depend on Rij, or equivalently, on Ru and 
Rjj separately. Fortunately, the dependence on Rjj is 
fully governed by Slater-Koster transformation rules^. 
Only one to three Slater-Koster integrals, depending on 
the angular momenta of ip^ and p v , arc needed to calcu- 
late the integral with any Rjj for fixed Rjj . These rules 
originate from the properties of spherical harmonics. 

The procedure is hence the following: we integrate nu- 
merically the required Slater-Koster integrals for a set of 
Ru, and store them in a table. This is done once for all 
orbital pairs. Then, for a given orbital pair, we interpo- 
late this table for Rjj , and use the Slater-Koster rules to 
get the overlap with any geometry — fast and accurately. 

Readers unfamiliar with the Slater-Koster transforma- 
tions can read the detailed discussion in Appendix [B] 
The numerical integration of the integrals is discussed in 
Appendix [Cj 

Before concluding this subsection, we make few re- 
marks about non-orthogonality. In DFTB it originates 
naturally and inevitably from Eq. (|48[) , because non- 
overlapping orbitals with diagonal overlap matrix would 
yield also diagonal Hamiltonian matrix, which would 
mean chemically non-interacting system. The transfer- 
ability of a tight-binding model is often attributed to 
non-orthogonality, because it accounts for the spatial na- 
ture of the orbitals more realistically. 

Non-orthogonality requires solving a generalized eigen- 
value problem, which is more demanding than normal 
eigenvalue problem. Non-orthogonality complicates, for 
instance, also gauge transformations, because the phase 
from the transformations is not well defined for the or- 
bitals due to overlap. The Pcicrls substitutional, while 
gauge invariant in orthogonal tight- bindin g 38 ' 39 , is not 



gauge invariant in non-orthogonal tight-binding (but af- 
fects only time-dependent formulation). 

C. Hamiltonian Matrix Elements 
From Eq. (|3~4")l the Hamiltonian matrix elements are 

K> = J (4 y2 + v M(r)) Mr), ( 49 ) 

where 

V s [n ] (r ) = V cxt {r) + V H [n ] (r) + V xc [n ] (r ) (50) 

is the effective potential evaluated at the (artificial) neu- 
tral density no(r) of the system. The density no is de- 
termined by the atoms in the system, and the above ma- 
trix element between basis states y, and v, in principle, 
depends on the positions of all atoms. However, since 
the integrand is a product of factors with three centers, 
two wave functions and one potential (and kinetic), all 
of which are non-zero in small spatial regions only, rea- 
sonable approximations can be made. 

First, for diagonal elements one can make a one- 
center approximation where the effective potential within 
volume Vi is 

V s [no](r) « V Si i[n ,i](r), (51) 

where /i G /. This integral is approximately equal to 
the eigenencrgics of free atom orbitals. This is only 
approximately correct since the orbitals <p^ are from the 
confined atom, but is a reasonable approximation that 
ensures the correct limit for free atoms. 

Second, for off-diagonal elements we make the two- 
center approximation: if fj, is localized around atom / and 
v is localized around atom J, the integrand is large when 
the potential is localized cither around / or J as well; 
we assume that the crystal field contribution from other 
atoms, when the integrand has three different localized 
centers, is small. Using this approximation the effective 
potential within volume Vi + Vj becomes 

V s [n }(r) « V sJ [n j](r) + V sJ [n ,j}(r), (52) 

where V s j[noj](r) is the Kohn-Sham potential with the 
density of a neutral atom. The Hamiltonian matrix ele- 
ment is 

H%, = J Mr)* (4v 2 + V sJ (r) + V a ,j(r)j Mr), 

(53) 

where fi £ I and v G J. Prior to calculating the integral, 
we have to apply the Hamiltonian to ip v . But in other 
respects the calculation is similar to overlap matrix ele- 
ments: Slater-Koster transformations apply, and only a 
few integrals have to be calculated numerically for each 
pair of orbitals, and stored in tables for future reference. 
See Appendix IC 21 for details of numerical integration of 
the Hamiltonian matrix elements. 
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IV. FITTING THE REPULSIVE POTENTIAL 

In this section we present a systematic approach to 
fit the repulsive functions V Tep (Ru) that appear in 
Eq. (|55|) — and a systematic way to describe the fitting. 
But first we discuss some difficulties related to the fitting 
process. 

The first and straightforward way of fitting is simple: 
calculate dimer curve Edft(R) for the element pair with 
DFT, require Edft(R) = Edftb(R), and solve 

V Icp (R) = E DFT (R) - [Ebs(R) + E coul (R)}. (54) 

We could use also other symmetric structures with N 
bonds having equal Ru ~ R, and require 

N ■ V Iep (R) = E DFT (R) - [E B s(R) + E coul (R)}. (55) 

In practice, unfortunately, it does not work out. The ap- 
proximations made in DFTB arc too crude, and hence 
a single system is insufficient to provide a robust repul- 
sion. As a result, fitting repulsive potentials is difficult 
task, and forms the most laborous part of parametrizing 
in DFTB. 

Let us compare things with DFT. As mentioned earlier, 
we tried to dump most of the difficult physics into the re- 
pulsive potential, and hence V rcp in DFTB has practical 
similarity to E xc in DFT. In DFTB, however, we have 
to make a new repulsion for each pair of atoms, so the 
testing and fitting labor compared to DFT functionals is 
multifold. Because xc-functionals in DFT are well doc- 
umented, DFT calculations of a reasonably documented 
article can be reproduced, whereas reproducing DFTB 
calculations is usually harder. Even if the repulsive func- 
tions are published, it would be a great advantage to be 
able to precisely describe the fitting process; repulsions 
could be more easily improved upon. 

Our starting point is a set of DFT structures, with 
geometries R, energies ^dft(-R), and forces F DFT (zero 
for optimized structures). A natural approach would be 
to fit V lcp so that energies ^dftb(-R) and forces F DFTB 
are as close to DFT ones as possible. In other words, 
we want to minimize force differences |i< lDFT — jr'DFTBj 
and energy differences I-Edft — -EdftbI on average for 
the set of structures. There are also other properties 
such as basis set quality (large overlap with DFT and 
DFTB wave functions), energy spectrum (similarity of 
DFT and DFTB density of states), or charge transfer to 
be compared with DFT, but these originate already from 
the electronic part, and should be modified by adjusting 
^confM and Hubbard U. Repulsion fitting is always the 
last step in the parametrizing, and affects only energies 
and forces. 

In practice we shall minimize, however, only force diffe- 
rences — we fit repulsion derivative, not repulsion directly. 
The fitting parameters, introduced shortly, can be ad- 
justed to get energy differences qualitatively right, but 
only forces are used in the practical fitting algorithm. 



There are several reasons for this. First, forces are ab- 
solute, energies only relative. For instance, since we do 
not consider spin, it is ambiguous whether to fit to DFT 
dimer curve with spin-polarized or spin-paired free atom 
energies. We could think that lower-level spin-paired 
DFT is the best DFTB can do, so we compare to spin- 
paired dimer curve — but we should fit to energetics of 
nature, not energetics in some flavors of DFT. Second, 
for faithful dynamics it is necessary to have right forces 
and right geometries of local energy minima; it is more 
important first to get local properties right, and after- 
wards look how the global properties, such as energy or- 
dering of different structural motifs, come out. Third, the 
energy in DFTB comes mostly from the band-structure 
part, not repulsion. This means that if already the band- 
structure part describes energy wrong, the short-ranged 
repulsions cannot make things right. For instance, if 
Euft(R) and £ , dftb(-R) for dimer deviates already with 
large R, short-range repulsion cannot cure the energetics 
anymore. For transferability repulsion has to be mono- 
tonic and smooth, and if repulsion is adjusted too rapidly 
catch up with DFT energetics, the forces will go wrong. 

For the set of DFT structures, we will hence minimize 
DFT and DFTB force differences, using the recipes be- 
low. 



A. Collecting Data 

To fit the derivative of the repulsion for element pair 
AB, we need a set of data points {Ri, V r ' ep (Ri)}. As men- 
tioned before, fitting to dimer curve alone does not give 
a robust repulsion, because the same curve is supposed 
to work in different chemical environments. Therefore it 
is necessary to collect the data points from several struc- 
tures, to get a representative average over different types 
of chemical bonds. Here we present examples on how to 
acquire data points. 



1. Force Curves and Equilibrium Systems 

This method can be applied to any system where all 
the bond lengths between the elements equal Rab or 
otherwise are beyond the selected cutoff radius i? cu t- In 
other words, the only energy component missing from 
these systems is the repulsion from N bonds between 
elements A and B with matching bond lengths. Hence, 



Ebftb{Rab) — Ebs{Rab) + E cou i(Rab) 
+E iep + N -V Iep (RAB), 



(56) 



where E vop is the repulsive energy independent of Rab ■ 
This setup allows us to change Rab, and we will require 



Vr ep (RAB) 



E^ ft (R A b) - [E' bs (R A b) + E'coui( r ab)} 



N 



(57) 



9 



where the prime stands for a derivative with respect to 
Rab- The easiest way is first to calculate the energy 
curve and use finite differences for derivatives. In fact, 
systems treated this way can have even different Rab's 
if only the ones that are equal are chosen to vary (e.g. 
a complex system with one appropriate AB bond on its 
surface). For each system, this gives a family of data 
points for the fitting; the number of points in the family 
does not affect fitting, as explained later. The dimer 
curve, with N = 1, is clearly one system where this 
method can be applied. For any equilibrium DFT struc- 
ture things simplify into 

t// /no \ ~E' BS (R AB ) — E' cou} (R AB ) 

V rep{ K AB) = Jr ! ( 58 i 

where R AB is the distance for which E^ FT (R AB ) = 0. 

2. Homonuclear Systems 

If a cluster or a solid has different bond lengths, the 
energy curve method above cannot be applied (unless 
a subset of bonds are selected). But if the system is 
homonuclear, the data points can be obtained the fol- 
lowing way. The force on atom / is 

Fj =F° + V; cp (Rij)Ru (59) 
=F°j + zijRiJ, (60) 

where F® is the force without repulsions. Then we min- 
imize the sum 

^I^DFTJ-F/I 2 (61) 
I 

with respect to e/j, with ejj = for pair distances larger 
than the cutoff. The minimization gives optimum e/j, 
which can be used directly, together with their Ri/s, as 
another family of data points in the fitting. 

3. Other Algorithms 

Fitting algorithms like the ones above are easy to con- 
struct, but a few general guidelines are good to keep in 
mind. 

While pseudo-atomic orbitals are calculated with 
LDA-DFT, the systems to fit the repulsive potential 
should be state-of-the-art calculations; all structural 
tendencies — whether right or wrong — are directly inher- 
ited by DFTB. Even reliable experimental structures can 
be used as fitting structures; there is no need to think 
DFTB should be parametrized only from theory. DFTB 
will not become any less density-functional by doing so. 

As data points are calculated by stretching selected 
bonds (or calculating static forces), also other bonds 



may stretch (dimer is one exception). These other bonds 
should be large enough to exclude repulsive interactions; 
otherwise fitting a repulsion between two elements may 
depend on repulsion between some other element pairs. 
While this is not illegal, the fitting process easily be- 
comes complicated. Sometimes the stretching can affect 
chemical interactions between elements not involved in 
the fitting; this is worth avoiding, but sometimes it may 
be inevitable. 



B. Fitting the Repulsive Potential 

Transferability requires the repulsion to be short- 
ranged, and we choose a cutoff radius i? cu t for which 
V lcp (R cut ) = 0, and also V r ' cp (R cut ) = for continuous 
forces. i?cut is one of the main parameters in the fitting 
process. Then, with given i? cu t, after having collected 
enough data points {Ri, V' we can fit the function 
V r ' ep (R). The repulsion itself is 

V Iep (R) = - V r ' cp (r)dr. (62) 

Jr 

Fitting of VJ.' ep using the recipe below provides a robust 
and unbiased fit to the given set of points, and the pro- 
cess is easy to control. We choose a standard smoothing 
spline^ for V^ ep (R) = U(R), i.e. we minimize the func- 
tional 

S [U(R)\ = J2 ( Kcp \ t/(J? ' ) ) 2 + A J Rmt U"(R) 2 dR 

(63) 

for total M data points {Ri, V T ' ep where U(R) is given 
by a cubic spline. Spline gives an unbiased representation 
for U (R) , and the smoothness can be directly controlled 
by the parameter A. Large A means expensive curva- 
ture and results in linear U(R) (quadratic y rop ) going 
through the data points only approximately, while small 
A considers curvature cheap and may result in a wiggled 
U(R) passing through the data points exactly. The pa- 
rameter A is the second parameter in the fitting process. 
Other choices for U(R) can be used, such as low-order 
polynomials^, but they sometimes behave surprisingly 
while continuously tuning i? cu t- For transferability the 
behavior of the derivative should be as smooth as possi- 
ble, preferably also monotonous (the example in Fig.[3Jt is 
slightly non-monotonous and should be improved upon). 

The parameters <ii are the data point uncertainties, 
and can be used to weight systems differently. With the 
dimension of force, oVs have also an intuitive meaning 
as force uncertainties, the lengths of force error bars. As 
described above, each system may produce a family of 
data points. We would like, however, the fitting to be 
independent of the number of points in each family; a fit 
with dimer force curve should yield the same result with 
10 or 100 points in the curve. Hence for each system 

Gi=a s ^fW s , (64) 
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V. ADJUSTING PARAMETERS 

In this section we summarize the parameters, give 
practical instructions for their adjustment, and give a 
demonstration of their usage. The purpose is to give an 
overall picture of the selection of knobs to turn while 
adjusting parametrization. 

Occasionally one finds published comments about the 
performance of DFTB. While DFTB shares flaws and 
failures characteristic to the method itself, it should be 
noted that DFTB paramctrizations are even more diverse 
than DFT functionals. A website in Ref. |U maintained 
by the original developers of the method, contains var- 
ious sets of paramctrizations. While these paramctriza- 
tions are of good quality, they are not unique. Namely, 
there exists no automated way of parametrizing, so that 
a straightforward process would give all parameters def- 
inite values. This is not necessarily a bad thing, since 
some handwork in parametrizing also gives feeling what 
is to be expected in the future, how well parametrizations 
are expected to perform. 



Pseudo- Atoms 



FIG. 3: (color online) a) Fitting the derivative of repulsive 
potential. Families of points from various structures, ob- 
tained by stretching C-H bonds and using Eq. (|57fl . Here 
Rcut — 1.8 A, for other details see Section [V] b) The repul- 
sive potential V lep (R), which is obtained from by integration 
of the curve in (a). 



where a s is the uncertainty given for system s, with N s 
points in the family. This means that systems with the 
same <7 s 's have the same significance in the process, ir- 
rcgardlcss of the number of data points in each system. 
The effect of Eq. is the same as putting a weight 
1/N S for each data point in the family. Note that A has 
nothing to do with the number of data points, and is 
more universal parameter. The cutoff is set by adding a 
data point at U(R cu t) = with a tiny a. 

Fig. [3] shows an example of fitting carbon-hydrogen 
repulsion. The parameters i? C ut and A, as well as param- 
eters <7j; , are in practice chosen to yield visually satisfying 
fit; the way of fitting should not affect the final result, and 
in this sense it is just a technical necessity — the simple 
objective is to get a smooth curve going nicely through 
the data points. Visualization of the data points can 
be also generally invaluable: deviation from a smooth 
behavior tells how well you should expect DFTB to per- 
form. If data points lay nicely along one curve, DFTB 
performs probably well, but scattered data points sug- 
gest, for instance, the need for improvements in the elec- 
tronic part. In the next section we discuss parameter 
adjustment further. 



The basis functions, and, consequently, the matrix 
elements arc determined by the confinement potentials 
14onf(0 containing the parameters ro for each element. 
The value ro = 2 • r cov , where r cov is the covalent radius, 
can be used as a rule of thumrA Since the covalent ra- 
dius is a measure for binding range, it is plausible that 
the range for environmental confining potential should 
depend on this scale — the number 2 in this rule is empir- 
ical. 

With this rule of thumb as a starting point, the quality 
of the basis functions can be inspected and ro adjusted by 
looking at (i) band-structure (for solids), (ii) densities of 
states (dimer or other simple molecules), or (iii) amount 
of data point scatter in repulsion fit (see Subsection llV B[ 
especially Fig. [3]) . Systems with charge transfer should 
be avoided herein, since the properties listed above would 
depend on electrostatics as well, which complicates the 
process; adjustment of ro should be independent of elec- 
trostatics. 

The inspections above are easiest to make with 
homonuclear interactions, even though hctcronuclear in- 
teractions are more important for some elements, such as 
for hydrogen. Different chemical environments can affect 
the optimum value of ro, but usually it is fixed for all 
interactions of a given element. 

B. Electrostatics 

Electrostatic energetics, as described in Subsec- 
tion HlDl are determined by the Hubbard U parameter, 
having the default value U = IE — EA. Since U is a 
value for a free atom, while atoms in molecules are not 
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free, it is permissible to adjust U in order to improve 
(i) charge transfer, (ii) density of states, (hi) molecular 
ionization energies and electron affinities, or (iv) excita- 
tion spectra for selected systems. Since U is an atomic 
property, one should beware when using several differ- 
ent elements in fitting — charge transfer depends on Z7's 
of other elements as well. 

Eq. relates U and FWHM of given atom together. 
But since Eq. (|2ip contains, in principle, also xc-contribu- 
tions, the relation can be relaxed, if necessary. Since 
FWHM affects only pair-interactions, it is better to ad- 
just the interactions directly like 

Cu -> Cu/xu, (65) 

where xu (being close to one) effectively scales both 
FWHM/ and FWHMj. If atomic FWHMs would be 
changed directly, it would affect all interactions and com- 
plicate the adjusting process. Note that FWHM affects 
only nearest neighbor interactions (see Fig. [T|) and al- 
ready next-nearest neighbors have (very closely) the pure 
1/i? interaction. 

An important principle, general for all parameters but 
particularly for electrostatics, is this: all parameters 
should be adjusted within reasonable limits. This means 
that, since all parameters have a physical meaning, if a 
parameter is adjusted beyond a reasonable and physically 
motivated limit, the parametrization will in general not 
be transferable. If a good fit should require overly large 
parameter adjustments (precise ranges are hard to give), 
the original problem probably lies in the foundations of 
tight-binding. 

C. Repulsive Potentials 

The last step in the parametrization is to fit the repul- 
sive potential. Any post-adjustment of other parameters 
calls for re-fitting of the repulsive potential. 

The most decisive part in the fitting is choosing the set 
of structures and bonds to fit. Parameters R cu t, <J and 
A are necessary, but they play only a limited part in the 
quality of the fit — the quality and transferability is de- 
termined by band structure energy, electrostatic energy, 
and the chosen structures. In fact, the repulsion fitting 
was designed such that the user has as little space for 
adjustment as possible. 

The set of structures should contain the fitted inter- 
action in different circumstances, with (i) different bond 
lengths, (ii) varying coordination, and (iii) varying charge 
transfer. In particular, if charge transfer is important for 
the systems of interest, calculation of charged molecules 
is recommended. 

A reasonable initial guess for the cutoff radius is i? cu t = 
1.5 x -Rdimcr, being half-way between nearest and next- 
nearest neighbors for homonuclear systems. It is then 
adjusted to yield a satisfying fit for the derivative of the 
repulsion, as in Fig. [31 while remembering that it has to 
be short ranged (i? cu t = 2 x i?dimer, for instance, is too 



large, lacks physical motivation, and makes fitting hard). 
Increasing R cut will increase V rep (R) at given R < i? cu t, 
which is an aspect that can be used to adjust energies 
(but not much forces). Usually the parameters a are used 
in the sense of relative weights between systems, as often 
they cannot be determined in the sense of absolute force 
uncertainties. The absolute values do not even matter, 
since the scale of a's merely sets the scale for the smooth- 
ness parameter A (you can start with a = 1 for the first 
system; if you multiply ct's by x, the same fit is obtained 
with A multiplied by x 2 ) — this is why the parameter A is 
not given any guidelines here. Large ct's can be used to 
give less weight for systems with (i) marginal importance 
for systems of interest, (ii) inter-dependence on other pa- 
rameters (dependence on other repulsions, on electrostat- 
ics, or on other chemical interactions), (iii) statistically 
peculiar sticking out from the other systems (reflecting 
situation that cannot be described by tight-binding or 
the urge to improve electronic part). 

D. Hydrocarbon Parametrization 

To demonstrate the usage of the parameters, we 
present the hydrocarbon parametrization used in this ar- 
ticle (this was first shot fitting without extensive adjust- 
ment, but works reasonably well). The Hubbard U is 
given by Eq. |17|) and FWHM by Eq. (J29J), both for hy- 
drogen and carbon. The force curves have been calcu- 
lated with GPAWi^ using the PBE xc-functionalii. 

Hydrogen: r = 1.08, U = 0.395. Carbon: r = 2.67, 
U = 0.376. Hydrogen-carbon repulsion: R cut = 3.40, 
A = 35, and systems with force curve: CH~, ethync, 
methane, benzene; all with ct; = 1. Carbon-carbon re- 
pulsion: R cut = 3.80, A = 200, and systems with force 
curve: C 2 , CC 2 ~, C 3 , C\~ with a t = 1, and C|~ with 
<Ti = 0.3. 

VI. EXTERNAL FIELDS 

Including external potentials to the formalism is 
straightforward: just add one more term in the Hamilto- 
nian of Eq. (|53|) . Matrix element with external (scalar) 
potential becomes, with plausible approximations, 

> + J ^(r)V ext (r)(p v (r)d 3 r 
w H^ v + V cxt (Ri) / ip*(p v + V ext {Rj) / 

* + \ {VL + V c J xt ) (66) 

for a smoothly varying V cxt (r). The electrostatic part in 
the Hamiltonian is 

and naturally extends Eq. (|4"2"1) . 



12 



VII. VAN DER WAALS FORCES 

Accurate DFT xc-functionals, which automatically 
yield the i?~ 6 long range attractive van der Waals inter- 
actions, are notoriously hard to make^. Since DFT in 
other respects is accurate with short-range interactions, 
it would be wrong to add van der Waals interactions by 
hand — addition inevitably modifies short-range parts as 
well. 

DFTB, on the other hand, is more approximate, and 
adding physically motivated terms by hand is easier. In 
fact, van der Waals forces in DFTB can conceptually be 
thought of as modifications of the repulsive potential. 
Since dispersion forces are due to a;c-contributions, one 
can see that for neutral systems, where 6n(r) = 0, disper- 
sion has to come from Eq. (jTUJ) . However, in practice it is 



better to leave V/ e p's short-ranged and add the dispersive 
forces as additional terms 



E 



vdW — 



E 

KJ 



(jlj 

flj(Rl.j)-^6~ 
U IJ 



(68) 



in the total energy expression. Here f(R) is a damping 
function with the properties 



f(R) = 



l, R>Ro 
0, R < R , 



(69) 



because the idea is to switch off van der Waals inter- 
actions for distances smaller than Rq, a characteristic 
distance where chemical interactions begin to emerge. 

The Cg-parameters depend mainly on atomic polariz- 
abilitics and have nothing to do with DFTB formalism. 
Care is required to avoid large repulsive forces, coming 
from abrupt behavior in f(R) near R ss Rq, which could 
result in local energy minima. For a detailed descriptions 
about the C6-parameters and the form of f(R) we refer 
to original Refs. and [U in this section we merely 
demonstrate how straightforward it is, in principle, to 
include van der Waals forces in DFTB. 



where u a (k,r) is function with crystal periodicity^ 9 . 
This means that a wave function ip a (k,r) changes by a 

phase e in translation T. We define new basis func- 
tions, not as localized orbitals anymore, but as Bloch 
waves extended throughout the whole crystal 



(fc,r) = ^£e* T ^(r-T), (71) 



where N is the (infinite) number of unit cells in the crys- 
tal. The eigenfunction ansatz 



(72) 



is then also an extended Bloch wave, as required by Bloch 
theorem, because k is the same for all basis states. Ma- 
trix elements in this new basis are 



and 



5^(fe,fe') = 5(k ~ k')S^(k) 



(73) 
(74) 



where 



S MU (k) = £ e lk - T (j V * (r)^(r - T) 
T 



(75) 



and similarly for H. Obviously the Hamiltonian con- 
serves k — that is why k labels the eigenstates in the first 
place. Note that the new basis functions are usually not 
normalized. 

Inserting the trial wave function (|72|) into Eq. (|38|) and 
by using the variational principle we obtain the secular 
equation 

J2 c a u (k) [JJ„„(fc) - e a (k)S^(k)} = 0, (76) 



VIII. PERIODIC BOUNDARY CONDITIONS 
A. Bravais Lattices 

Calculation of isolated molecules with DFTB is 
straightforward, but implementation of periodic bound- 
ary conditions and calculation of electronic band- 
structures is also easy^. As mentioned in the introduc- 
tion, this is usually the first encounter with tight-binding 
models for most physicists; our choice was to discuss pe- 
riodic systems at later stage. 

In a crystal periodic in translations T, the wave func- 
tions have the Bloch form 



1pa(k, 



e lk - r u a {k,r) 



(70) 



where 



H llv {k)=Hl v {k) + h} llv S llv {k) 



(77) 



for each fc-point from a chosen set, such as Monkhorst- 
Pack sampled^. Above we have 



1 



(78) 



as in Eq. |42|) . and Mullikcn charges arc extensions of 
Eq. (E) , 



^ = EE^( fc ) E hc a ;(kK(k)s^(k) +c.c.} 



a k 



(79) 
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The sum for the electrostatic energy per unit cell, 

^ unit cell 

E coul = - E Y^1ij(Rij-T)A qi Aqj, (80) 
u T 

can be calculated with standard methods, such as Ewald 
summational, and the repulsive part, 

unit cell 

J2Vrlt(Ru) = 2 E Y, v "p( R "~ T ) ( 81 ) 

KJ IJ T 

is easy because repulsions are short-ranged (and 
V le p(0) = is understood). 

B. General Symmetries 

Thanks to the transparent formalism of DFTB, it is 
easy to construct more flexible boundary conditions, such 
as the "wedge boundary condition" introduced in Ref. 
This is one example of DFTB in method development. 

General triclinic unit cells arc copied by translations, 
and DFT implementation is easy with plane waves, real- 
space grids or localized orbitals with fixed quantization 
axis. But if we allow the quantization axis of localized 
orbitals to be position-dependent, we can treat more 
general symmetries which have rotational symmetry^ 
or even combined rotational and translational (chiral) 
symmetries 5 ^. 

The basic idea is to enforce the orbitals to have the the 
same symmetry as the system. This requires that basis 
functions not only depend on atom positions like 

Mr -A,,), (82) 
as usual, but more generally like 

D(i? M )Mr), (83) 

where D{R lJ ) is an operator transforming the orbitals in 
any position-dependent manner, including both transla- 
tions and rotations. The only requirement is that the 
orbitals are complete and orthonormal for a given angu- 
lar momentum. If the quantization axes change, things 
become unfortunately messy. However, suitably defined 
basis orbitals yield well-defined Hamiltonian and over- 
lap matrices, and enable simulations of systems like bent 
tubes or slabs, helical structures such as DNA, or a 
piece of spherical surface — with a greatly reduced num- 
ber of atoms. Similar concepts are familiar from chem- 
istry, where symmetry-adapted molecular orbitals are 
constructed from the atomic orbitals, and computational 
effort is hereby reduced 55 . Detailed treatment of these 
general symmetries is a work in progress^. 

IX. DENSITY-MATRIX FORMULATION 

In this section we introduce DFTB using density- 
matrix formulation. We do this because not only does 



the formulation simplify expressions, but it also makes 
calculations faster in practice. This practical advantage 
comes because the density-matrix, 

/v = E /«(<£<$•) = E f°W> ( 84 ) 

a a 

contains a loop over eigenstates; quantities calculated 
with Pf^u simply avoid this extra loop. It has the proper- 
ties 

pSp = p (~ idempotency) (85) 

P t tv = P*u li (p = P t ) (86) 
Tr (p a S) = 1 (eigenfunction normalization). (87) 

Wc define also the energy-weighted density-matrix 

a 

and symmetrized density matrix 

Pi»> = \{P^ + M)> (89) 

which is symmetric and real. Using p^ we obtain simple 
expressions, for example, for 

E BS = Tr{pH°) (90) 

N el = Tr(pS) = E E Pv S "i» ( 91 ) 

Ql = Tr/(pS) = E E p^« s ^ ( 92 ) 

where Ebs is the band-structure energy, N e i is the total 
number of electrons, and qi is the Mulliken population 
on atom /. Tr/ is partial trace over orbitals of atom / 
alone. 

It is practical to define also matrices' gradients. They 
do not directly relate to density matrix formulation, but 
equally simplify notation, and are useful in practical im- 
plementations. We define (with Vj = d/dRj) 

dS^y = V,/5, liy p, G /, v e J 

= J Mr - H/)VjV„(r - Rj) (93) 

= (MVjM> 

and 

dH^ = VjH^ v pel, v e J 
= (m#|VjM 

with the properties dS^ v = —dS ufl and dH^ = 
dH*^. From these definitions we can calculate ana- 
lytically, for instance, the time derivative of the overlap 
matrix for a system in motion 

S=[dS,V], (95) 



14 



with commutator [A, B] and matrix V ^ v ~ S^Ri, 
fx £ I. Force from the band-energy part, the first line 
in Eq. (|43]l . for atom / can be expressed as 

Fj = -Tr I (pdH - p e dS) + c.c, (96) 

which is, besides compact, useful in implementation. The 
density-matrix formulation introduced here is particu- 
larly useful in electronic structure analysis, discussed in 
the following section. 

X. SIMPLISTIC ELECTRONIC STRUCTURE 
ANALYSIS 

One great benefit of tight-binding is the ease in ana- 
lyzing the electronic structure. In this section we present 
selected analysis tools, some old and renowned, others ca- 
sual but intuitive. Other simple tools for chemical anal- 
ysis of bonding can be found from Ref. IHtI 

A. Partial Mulliken Populations 

The Mulliken population on atom J, 

qi = Tr/(p5) = E E P^ S '^ = E^w ( 97 ) 
fj.ei v fj.ei 

is easy to partition into smaller pieces. Population of a 
single orbital fi is 

?M = E P^ Sv » = O 55 )^ ( 98 ) 

while population on atom I due to eigenstate i\) a alone is 

qi,a = Ti 7 (p Q 5) = E E %»> S »l» (") 

so that 

E (j2 /a«f,-l = E>) = N * ( 10 °) 

Population on orbitals of atom / with angular momentum 
I is, similarly, 

q\= E ( 101 ) 

The partial Mulliken populations introduced above are 
simple, but enable surprisingly rich analysis of the elec- 
tronic structure, as demonstrated below. 

B. Analysis Beyond Mulliken Charges 

At this point, after discussing Mulliken population 
analysis, we comment on the role of wave functions in 



DFTB. Namely, internally DFTB formalism uses atom 
resolution for any quantity, and the tight-binding spirit 
means that the matrix elements iT M „ and S M „ are just pa- 
rameters, nothing more. Nonetheless, the elements H^ u 
and Sfa, are obtained from genuine basis orbitals <p m (t) 
using well-defined procedure — these basis orbitals remain 
constantly available for deeper analysis. The wave func- 
tions are 

Mr) = J2 c >^ ( 102 ) 
and the total electron density is 

n(r) = E fa\Hr)\ 2 = E P^'MM*)* (103) 

awaiting for inspection with tools familiar from DFT. 
One should, however, use the wave functions only for 
analysis^. The formalism itself is better off with Mul- 
liken charges. But for visualization and for gaining 
understanding this is a useful possibility. This distin- 
guishes DFTB from semiempirical methods, which — in 
principle — do not possess wave functions but only ma- 
trix elements (unless made ad hoc by hand). 

C. Densities of States 

Mulliken populations provide intuitive tools to inspect 
electronic structure. Let us first break down the energy 
spectrum into various components. The complete energy 
spectrum is given by the density of states (DOS), 

DOS(e) = E>(e-ea), (104) 

a 

where S a (e) can be either the peaked Dirac delta- 
function, or some function — such as a Gaussian or a 
Lorentzian — with broadening parameter a . DOS carry- 
ing spatial information is the local density of states, 

LDOS(e,r) =E^(£"£a)|^a(r)| 2 , (105) 

a 

with integration over J d 3 r yielding DOS (e). Sometimes 
LDOS(r) = E/>a(r)| 2 , (106) 

a 

where f' a are weights chosen to select states with 
given energies, as in scanning tunneling microscopy 
simulations^. Mulliken charges, pertinent to DFTB, 
yield LDOS with atom resolution, 

LD0S(e,I)=£V(e-e o ) 9ZlO , (1° 7 ) 

a 

which can be used to project density for group of atoms 
TZ as 

LDOS^(e) = E LDOS(e, I). (108) 
ien 
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TABLE I: Simplistic electronic structure and bonding analysis 
for selected systems: C2-H2 (triple CC bond), C2-H4 (double 
CC bond), C2H6 (single CC bond), benzene, and graphene 
(r-point calculation with 64 atoms). We list most energy and 
bonding measures introduced in the main text (energies in 
eV). 



property 


C2H2 


C2H4 


C2-H6 


benzene 


graphene 


q A H 


0.85 


0.94 


0.96 


0.95 




A H 


1.23 


0.45 


0.26 


0.32 




ABh 


-2.75 


-3.27 


-3.34 


-3.39 




-Eprom,ff 


0.97 


0.41 


0.25 


0.30 






4.15 


4.13 


4.12 


4.05 


4.00 


A c 


5.86 


5.86 


5.78 


6.48 


6.94 


AB C 


-9.16 


-9.74 


-10.45 


-9.74 


-9.78 


-Eprom,C 


5.62 


5.69 


5.64 


6.46 


6.94 


B C H 


-8.14 


-7.90 


-7.84 


-7.93 




Bcc 


-22.01 


-15.82 


-9.39 


-13.01 


-12.16 


M C H 


0.96 


0.95 


0.97 


0.96 




Mac 


2.96 


2.02 


1.01 


1.42 


1.25 



For instance, if systems consists of surface and adsorbed 
molecule, we can plot LDOS mo i(e) and LDOS sur f(e) to 
see how states are distributed; naturally LDOS mo i(e) + 
LDOS surf (e) =DOS(e). 

Similar recipes apply for projected density of states, 
where DOS is broken into angular momentum compo- 
nents, 

PDOS(£,0=E<H £ -^)E^a' ( 109 ) 
a J 

such that, again J2i PDOS(e, I) = DOS(e). 



D. Mayer Bond-Order 

Bond strengths between atoms are invaluable chemi- 
cal information. Bond order is a dimcnsionlcss number 
attached to the bond between two atoms, counting the 
differences of electron pairs on bonding and antibonding 
orbitals; ideally it is one for single, two for double, and 
three for triple bonds. In principle, any bond strength 
measure is equally arbitrary; in practice, some measures 
are better than others. A measure suitable for many pur- 
poses in DFTB is Mayer bond-order—, defined for bond 
I J as 

Mu= ]T (pS)^(pS)^. (110) 

The off-diagonal elements of pS can be understood as 
Mulliken overlap populations, counting the number of 
electrons in the overlap region — the bonding region. It 
is straightforward, if necessary, to partition Eq. (jllOp 
into pieces, for inspecting angular momenta or eigenstate 
contributions in bonding. Look at Refs. [6(3 and for 
further details, and Table U for examples of usage. 



E. Covalent Bond Energy 

Another useful bonding measure is the covalent bond 
energy, which is not just a dimcnsionlcss number but 
measures bonding directly using energy^. 

Let us start by discussing promotion energy. When free 
atoms coalesce to form molecules, higher energy orbitals 
get occupied — electrons get promoted to higher orbitals. 
This leads to the natural definition 

Eprom = y^(g( M ) ~ 9(™) 6 at ° m )#° M - (HI) 

Promotion energy is the price atoms have to pay to pre- 
pare themselves for bonding. Noble gas atoms, for in- 
stance, cannot bind to other atoms, because the promo- 
tion energy is too high due to the large energy gap; any 
noble gas atom could in principle promote electrons to 
closest s-statc, but the gain from bonding compared to 
the cost in promotion is too small. 

Covalent bond energy, on the other hand, is the energy 
reduction from bonding. We define covalent bond energy 
as£i 

E C ov = {EbS — Ef ICC atoms) _ E prom . (H2) 

This definition can be understood as follows. The term 
(Ebs — E[ ICC atoms) is the total gain in band-structure 
energy as atoms coalesce; but atoms themselves have to 
pay -Eprom, an on-site price that does not enhance binding 
itself. Subtraction gives the gain the system gets in bond 
energies as it binds together. More explicitly, 

E cov = E B s - ^ Q{n)H ^ (113) 
= J2p flv (H°„-e^S> lv ), (114) 

where 

This can be resolved with respect to orbital pairs and 
energy as 

E cov ^(e) = ]T 6°(e - e a )p%{H^ - e Vfl S Vfi ). (116) 

a 

E C av,n,v can be viewed as the bond strength between or- 
bitals p, and v — with strength directly measured in en- 
ergy; negative energy means bonding and positive energy 
antibonding contributions. Sum over atom orbitals yields 
bond strength information for atom pairs 

E covJ j(e) = ^^(-E COV)AU/ (e) + c.c), (117) 
and sum over angular momentum pairs 

l»=la l u =h 

E l ±{e) = J2 E (£cov,^(e) + [cc.]) (118) 

fl V 
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FIG. 4: (color online) Covalent bonding energy contributions 
in graphene (F-point calculation with 64 atoms in unit cell). 
Both bonding and antibonding ss bonds are occupied, but 
sp and pp have only bonding contributions. At Fermi-level 
(zero-energy) the pp-bonding (the 7r-cloud above and below 
graphene) transforms into antibonding (n*) — hence any ad- 
dition or removal of electrons weakens the bonds. 



gives bonding between states with angular momenta l a 
and lb (no c.c. for l a = It). For illustration, we plot 
covalent bonding contributions in graphene in Fig. [4] 



and 

BlJ =V x ll + 1IJ Aq I Aqj 

+ J2 p^( H °u~ s ^)+ c - c - (124) 

We call Ai the absolute atom energy of atom J, and Bjj 
the absolute bond energy of bond IJ. Measuring atom 
energies with Aj and bond energies with Bjj explicitly 
takes into account all energetics. Eq. ()121|) can be further 
simplified into 

e ' = Y,\ a i + \Y. b A =J2 ab ^ ( 125 ) 

where ABj measures how much atom I contributes to 
total binding energy — in electrostatic, repulsive, promo- 
tive, and bonding sense. For homonuclear crystals the 
binding energy per atom is directly ABj, and for het- 
eronuclear systems the binding energy per atom (the 
negative of cohesion energy) is averaged ABj; positive 
ABj means atom / would rather be a free atom, even 
though for charged systems the interpretation of these 
numbers is more complicated. Visualizing Ai, Bu, and 
AB\ gives a thorough and intuitive measure of energetics; 
see Table U for illustrative examples. Note that Aj, Bjj, 
and ABi come naturally from the exact DFTB energy 
expression — they are not arbitrary definitions. 



F. Absolute Atom and Bond Energies 

While Mayer bond-order and E cov are general tools for 
any tight-binding flavor, neither of them take electrostat- 
ics or repulsions between atoms into account. Hence, to 
conclude this section, we introduce a new analysis tool 
that takes also these contributions into account. 

The DFTB energy with subtracted free-atom energies, 

E =E — Ef ICC atoms (119) 

=Tr(pH°) + ^ 7 „A^A^ (120) 
ij 

i \ * \rU \ * „frcc atom tj 
KJ n 

can be rearranged as 

E' = Y j Ai + Y j Bi,j, (121) 

/ KJ 

where 

Aj = i 7// A 9/ 2 + - qf;f atom )^ (122) 

= pn&qi + E piomJ (123) 



XI. CONCLUSIONS 

Here ends our journey with DFTB for now. The road 
up to this point may have been long, but the contents 
have made it worth the effort: we have given a detailed 
description of a method to do realistic electronic struc- 
ture calculations. Especially the transparent chemistry 
and ease of analysis makes DFTB an appealing method 
to support DFT simulations. With these features tight- 
binding methods will certainly remain an invaluable sup- 
porting method for years to come. 
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APPENDIX A: CALCULATING THE DFT 
PSEUDO-ATOM 

The pseudo-atom, and also the free atom without the 
confinement, is calculated with LDA-DFT 6? . More re- 
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cent xc-functionals could be used, but they do not im- 
prove DFTB parametrizations, whereas LDA provides a 
fixed level of theory to build foundation. However, bet- 
ter DFT functionals can — and should be — used in the 
repulsive potential fitting; see Section ITVl 

Elements with small atomic numbers are calculated us- 
ing non-relativistic radial Schrodingcr equation. But for 
some elements, such as gold, chemistry is greatly modi- 
fied by relativistic effects, which have to be included in 
the atom calculation. 

In the four-component Dirac equation with central po- 
tential good quantum numbers are energy, total angular 
momentum j, its z-component j z , and — k which is the 
eigenvalue of the operator 



K = 



Si 





1 
-S • L - 1 



(Al) 



where L is the orbital angular momentum operator and 
the components of 4 x 4 relativistic spin-matrix S are 



_ ( <Jk 
a k 



(A2) 



with 2x2 Pauli spin-matrices <ik ■ Remember that angu- 
lar momentum I is not a good quantum number; the up- 
per and lower two components are separately eigenstates 
of L 2 with different angular momenta, and coupled by 
spin-orbit interaction. In other words, a given I (that we 
are interested in) appears in 4-componcnt spinors states 
with different j. 

The intention is to include relativistic effects from the 
Dirac equation, but still use the familiar non-relativistic 
machinery. This can be achieved by ignoring the spin- 
orbit interaction, decoupling upper and lower compo- 
nents. By considering only the upper components as a 
non-relativistic limit, I becomes again a good quantum 
number. The radial equation transforms into the scalar- 
relativistic equatio n 63 ' 64 



d 2 Rni( r) (1(1 + 1) 
dr 2 
1 



2M(r)(V s (r) - e nl ) R nl (r) 



dM(r) (dR nl (r) 



M (r) dr 



Rni{r) 



dr 



= 0. 



(A3) 



Here a = 1/137.036 is the fine structure constant, 

V.(r) = -| + V H (r) + V^ A (r) + V coni (r), (A4) 
with or without the confinement, and we defined 



C wz, 





FIG. 5: (color online) Illustrating overlap integral calcula- 
tion, (a) Originally one p x orbital locates at origin, another 
p x orbital at R. (b) Coordinate system is rotated so that an- 
other orbital shifts to Rz; this causes the orbitals in the new 
coordinate system to become linear combinations of p x and 
Pz. Hence the overlap can be calculated as the sum of the 
so-called S'(pp7r)-integral in (c), and S(ppo")-integral in (d). 
The integrals in (e) and (f) are zero by symmetry. 



trick, and is, for instance, used routinely for generating 
DFT pseudo-potentials^. 

The pseudo-atom calculation, as described, yields the 
localized basis orbitals we use to calculate the matrix ele- 
ments. Our conventions for the real angular part F^(#, if) 
of the orbitals <fin(r) = R fJ ,(r)Y IJj (9,ip) are shown in Ta- 
bic |TTJ The sign of R^{r) is chosen, as usual, such that 
the first antinodc is positive. 



APPENDIX B: SLATER-KOSTER 
TRANSFORMATIONS 



Consider calculating the overlap integral 



M(r) = l + —[e nl -V s (r)}. 



(A5) 



The reminiscent of the lower two components of Dirac 
equation is (k), which is the quantum number k averaged 
over states with different j, using the degeneracy weights 
j(j + 1); a straightforward calculation gives (k) = —1. 

Summarizing shortly, for given I the potential in the 
radial equation is the weighted average of the poten- 
tials in full Dirac theory, with ignored spin-orbit inter- 
action. This scalar-relativistic treatment is a standard 



S xx = J Px (r)Px(r - R)d 3 r (Bl) 



for orbital p x at origin and another p x orbital at _R, as 
shown in Fig. [5k.. We rotate the coordinate system pas- 
sively clockwise, such that the orbital previously at R 
shifts to R = Rz in the new coordinate system. Both 
orbitals become linear combinations of p x and p z in the 
new coordinate system, p x — > p x sin a + p z cos a, and the 
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TABLE II: Spherical functions, normalized to unity ( J Y(8, cp) sin ddddip — 1) and obtained from spherical harmonics as 
Y oc Yi m ±Y* m . Note that angular momentum I remains a good quantum number for all states, but magnetic quantum number 
m remains a good quantum number only for s-, p z -, and d 3z 2_ r 2-states. 



visualization 



symbol 



Y(8,<P) 




'Y 



o 



Py( e >v) y^sinesin^ 

Pz( e -.<P) \fi; cose 
d 3z 2_ r 2(6,ip) y^(3 cos 2 0- 



d xy (6,ip) 

d yz {6,ip) 
d zx (8,ip) 



sin 2 8 cos 2o3 

Ion r 



sin 2 8 sin 2ip 



sin 28 cos 



(B2) 



overlap integral becomes a sum of four terms 
J Px(r)p x (r - Rz) ■ sin 2 a (Fig. Efc) 
+ y Pz(r)p z {r - Rz) ■ cos 2 a (Fig. [SJi) 
+ y Px(r)p z (r — RZ) ■ cos a sin a (Fig. [5^) 
+ y Pzi^Pxir — Rz) ■ cos a sin a (Fig. [5f). 



The last two integrals are zero by symmetry, but the first 
two terms can be written as 

S xx = x 2 S(ppa) + (1 - x 2 )S(ppir), (B3) 

where x = cos a is the direction cosine of R. For sim- 
plicity we assumed y = 0, but the equation above applies 
also for y (using hindsight we wrote (1 — x 2 ) instead 
of z 2 ). The integrals 



S(ppa) = / p z (r)p z (r — Rz)d 3 r 



(B4) 



are called Slater-Koster integrals and Eq. (|B3[) is called 
the Slater-Koster transformation rule for the given or- 
bital pair (orbitals may have different radial parts; the 
notation S fJil/ (ppa) stands for radial functions R^{r) and 
R u {r) in the basis functions fi and v). Similar reason- 
ing can be applied for other combinations of p-orbitals 



as well — they all reduce to Slater-Koster transformation 
rules involving S(ppa) and S{ppir) integrals alone. This 
means that only two integrals with a fixed R is needed 
for all overlaps between any two p-orbitals from a given 
element pair. 

Finally, it turns out that 10 Slater-Koster integrals, 
labeled dda, ddn, ddS, pda, pdir, ppcr, ppir, sda, spa, and 
ssa, are needed to transform all s-, p-, and d- matrix 
elements. The last symbol in the notation, er, tt, or i5, 
refers to the angular momentum around the symmetry 
axis, and is generalized from the atomic notation s, p, d 
for I = 0,1,2. 

Tabic Hill shows how to select the angular parts for cal- 
culating these 10 Slater-Koster integrals. The integrals 
can be obtained in many ways; here we used our setup 
with the first orbital at the origin and the second or- 
bital at Rz. This means that we set the direction cosines 
2 = 1 and x = y = in Table HVl and chose orbital pairs 
accordingly. 

Finally, Table IIVI shows the rest of the Slater-Koster 
transformations. The overlap Sy. v between orbitals ip^ at 

R is the sum 



J? M = and (p v at R v 



Su 



c t S^ v [t), 



(B5) 



for the pertinent Slater-Koster integrals t (at most 
three). The gradients of the matrix elements come di- 
rectly from the above expression by chain rule: Slater- 
Koster integrals S^^r) depend only on R and the coef- 
ficients c T only on R. 

For 9 orbitals (one s-, three p-, and five d-orbitals) 
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TABLE III: Calculation of the 10 Slater-Koster integrals r. The first orbital (with angular part ti) is at origin and the 
second orbital (with angular part T2) at Rz. 4> T (Qi,Q2) is the function resulting from azimuthal integration of the real spherical 
harmonics (of Table Ell <j> T (O 1 ,0 2 ) = / =0 dipY^ ip)Y T2 (6 2 , <f), and is used in Eq. (JC6j) . The images in the first row visualize 
the setup: orbital (o) is centered at origin, and orbital (x) is centered at Rz; shown are wave function isosurfaces where the 
sign is color-coded, red (light grey) is positive and blue (dark grey) negative. 



ddS d x 




t y ri (0i lV >) y T2 (9 2 , v ) <M0i>02) 

dda d 3z 2_ r 2 d 3z 2_ r 2 |(3cos 2 9i — l)(3cos 2 6 2 — 1) 
dd-K d zx d zx 



ppn p x 



sda s 



spa s 



d x 



if sin 2 61 sin 2 9 2 



^ cos #1 (3 cos 2 6 2 - 1) 

sin 81 sin 62 cos 9 2 

I cos 6\ cos 6 2 

§ si: 
\/Si 



d 3z 3_ r 2 ^(3cos 2 6» 2 - 1) 
p z ^ cos 62 



81 transformations are required, whereas only 29 are in 
Table IIV1 Transformations with an asterisk can be ma- 
nipulated to yield another 16 transformations and the 
remaining ones can be obtained by inversion, which ef- 
fectively changes the order of the orbitals. This inversion 
changes the sign of the integral according to orbitals' an- 
gular momenta and l v , 



S i1v {t) = S vim {t)-{-1) 1 » +1 % (B6) 



because orbital parity is (— 1) . 

The discussion here concentrated only on overlap, 
but the Slater-Koster transformations apply equally for 
Hamiltonian matrix elements. 



APPENDIX C: CALCULATING THE 
SLATER-KOSTER INTEGRALS 

1. Overlap Integrals 

We want to calculate the Slater-koster integral 

S» v (t) = {ippTilVur,}, (CI) 

with 

if^ (r) = R»{r)Y Tl (6, tp) = R^(r)Q Tl (0)$ T1 (<p) (C2) 
and 

VurM = Ru(r)Y n (e,<p) = R„(r)e T2 (6)$ T2 (ip), (C3) 

where R(r) is the radial function, and the angular parts 
Y Ti (9, if) are chosen from Table IIIII and depend on the 
Slater-Koster integral r in question. Like in Appendix IBl 
we choose ip^ to be at the origin, and <p v to be at R = Rz. 
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TABLE IV: Slater-Koster transformations for s-, p-, and d-orbitals, as first published in Ref. |36|. To shorten the notation 
we used a = x 2 + y 2 and (3 = x 2 — y 2 . Here x, y, and z are the direction cosines of R, with x 2 + y 2 + z 2 = 1. Missing 
transformations are obtained by manipulating transformations having an asterisk (*) in the third column, or by inversion. Here 
7Ji is <^ M 's angular momentum and n is y>„'s angular momentum. 



p at 


v at R 


s 


s 






Px 




s 






s 


d x 2 _ 


y 


s 


d 3z 2 


-r 2 


Px 


pi- 




Px 


Py 




Px 


Pz 




Px 


dxy 




Px 


d y 2 




Px 


dzx 




Px 


d x 2 - 


y 2 


i y 


d x 2_ 


y2 


Pz 


d x 2_ 


y 2 


Px 


d r iz 2 


-r 2 


i y 


d 3z 2 


-r 2 


p z 


d 3z 2 


-T 2 


dxy 








dyz 




dxy 


dzx 




dxy 


d x 2_ 


y 2 


dyz 


d x 2_ 


y 2 


dzx 


d x 2_ 


y 2 


dxy 


d-iz 2 


-r 2 


dyz 


d 3 z 2 


-r 2 


dzx 


d 3 z 2 


-r 2 


d x 2 —y 2 


dx 2 ~ 


I/ 2 


d x 2 — y 2 


d-iz 2 


-r 2 


d 3z 2_ r 2 


d 3 z 2 


-r 2 



~T 

x 

V3xy 

2 1 
Z — jQ 

x 2 1-x 2 

xy —xy 

xz —xz 

V3x 2 y 2/(1 - 2a; 2 ) 

\J?>xyz —2xyz 

V3x 2 z z(l - 2x 2 ) 

±V3xf3 x(l-/3) 

§V%? -VQ- + P) 



xz 2 



\\F&zfi -z/3 

x(z 2 — ^a) — v3: 

y(z 2 - \a) -V3yz 2 

z(z 2 — |a) \J~3zot 

3x 2 y 2 ' a - 4x 2 y 2 z 2 + x 2 y 2 

3xy 2 z xz(l — Ay 2 ) xz(y 2 — 1) 

3x 2 yz yz{l — 4x 2 ) yzix 2 — 1) 

lxy/3 -2xyf3 ±xy[3 

hzf3 -yz(l + 2p) V z{l + \P) 

\zxfi zx(l - 2/3) -xz(l - 

^/3xy(z 2 — |a) — 2\fixyz 2 ^y/3xy{l + z 2 

V3yz(z 2 — ~a) V3yz(a — z 2 ) —\V3yza 

\/3xz(z 2 — |a) -\/3~xz(a — z 2 ) ~^\/3xza 

-f3 2 a — f3 2 z 2 + -0 1 

tV3/3{z 2 - \a) -V3z 2 p + z 2 )/3 

(z 2 - ±a) 2 3z 2 a fa 2 



Explicitly, 

S^(t) = J ^d 3 r[^(ri)e Tl (ei)$ Tl (^i; 



(C4) 



x [R„{r 2 )e T2 (9 2 )$ T2 (tp 2 )}, 



where r\ = r and r 2 = r — Rz. Switching to cylindrical 
coordinates we get 



S^{t) = J J dzpdpi? A1 (ri)i?^(r 2 ) 



xe T1 (^)e T2 (0 2 ) / $ T1 (^ 1 )<i> T2 (v 32 )d^, 

(C5) 

and we see that since z-axis remains the symmetry axis, 
the (^-integration can be done analytically. The sec- 
ond line in Eq. (|C5[) becomes an analytical expression 
4>t {01,02), and is given in Table HTT1 Note here that r is 
a spherical coordinate, whereas p is the distance from a 
£-axis in cylindrical coordinates. We are left with 



dzpdpR ll {r 1 )R„{r 2 )M0i,02), (C6) 



a two-dimensional integral to be integrated numerically. 



2. Hamiltonian Integrals 

The calculation of the Slater-Koster Hamiltonian ma- 
trix elements 



(C7) 



is mostly similar to overlap. The potentials V s j[noj]{r) 
in the Hamiltonian 

H° = -\v 2 + V s j[n j](r) + V s ,j[n 0jJ ]{r), (C8) 
with p € I and v G J, arc approximated as 

K,iK/](r) « V'T^r) - V coniJ (r), (C9) 

where V^ 111 ^) is the self-consistent effective potential 
from the confined pseudo-atom, but without the confin- 
ing potential. The reasoning behind this is that while 
14onf(0 yields the pseudo-atom and the pseudo-atomic 
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orbitals, the potential V s [no](r) in H° should be the ap- 
proximation to the true crystal potential, and should not 
be augmented by confinements anymore. The Hamilto- 
nian becomes 



The form used in numerical integration is 



* s,J 



r) - V con { t j(r) 



(CIO) 



and the matrix element 

+ {<Pm \V s m f l - V con{ j - V conil j\<p„ T2 ) (Cll) 
=£y° S^ v (r) 

+ {PuT^Vgj 1 — Konf,I — Vconf,j\(PuT 2 } j (C12) 

depending whether we operate left with — V 2 /2 + VT? 

or right with — V 2 /2 + VT} nf ; ip^s are eigenstates of the 

confined atoms with eigenvalues £^ (including the con- 
finement energy contribution which is then subtracted). 



dzpdpR^ (r i ) R v (r 2 ) <\> T (#1 , 82 ) 



x [V s c T\n) - Konf,i(ri) - Konf,j(r 2 )] . 

(C13) 

As an internal consistency check, we can operate to 
tp„ also directly with V 2 , which in the end requires 
just d 2 R l/ {r)/dr 2 , but gives otherwise similar integration. 
Comparing the numerical results from this and the two 
versions of Eqs. (|C11[) and (|C12[) give way to estimate 
the accuracy of the numerical integration. 

Note that the potential in Eq. (jCllj) diverges as r — > 
Rz, and the potential in Eq. (|C11|) diverges as r — > 0. 
For this reason we use two-center polar grid, centered at 
the origin and at Rz, where the two grids are divided by 
a plane parallel to xy-planc, and intersecting with the 
z-axis at ^R ■ z. 



t Electronic address: pekka.koskinen@iki.fi 

1 A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. 
Novoselov, A. K. Geim, The electronic properties of 
graphene, Rev. Mod. Phys. 81 (2009) 109. 

2 D. Porezag, T. Frauenheim, T. Kohler, G. Seifert, 
R. Kaschner, Construction of tight-binding-like potentials 
on the basis of density-functional theory: application to 
carbon, Phys. Rev. B 51 (1995) 12947. 

3 M. Elstner, D. Porezag, G. Jungnickel, J. Elstner, 
M. Haugk, T. Frauenheim, S. Suhai, G. Seifert, Self- 
consistent-charge density-functional tight-binding method 
for simulations of complex materials properties, Phys. Rev. 
B 58 (1998) 7260. 

4 M. Elstner, T. Frauenheim, E. Kaxiras, G. Seifert, 
S. Suhai, A self-consistent charge density- functional based 
tight-binding scheme for large biomolecules, phys. stat. sol. 
(b) 217 (2000) 357. 

5 T. Frauenheim, G. Seifert, M. Elstner, Z. Hajnal, G. Jung- 
nickel, D. Porezag, S. Suhai, R. Scholz, A self-consistent 
charge density-functional based tight-binding method for 
predictive materials simulations in physics, chemistry and 
biology, phys. stat. sol. b 217 (2000) 41. 

6 P. Koskinen, H. Hakkinen, B. Huber, B. von Issendorff, 
M. Moseler, Liquid-liquid phase coexistence in gold clus- 
ters: 2d or not 2d?, Phys. Rev. Lett. 98 (2007) 015701. 

7 K. A. Jackson, M. Horoi, I. Chaudhuri, T. Frauenheim, 
A. A. Shvartsburg, Unraveling the shape transformation 
in silicon clusters, Phys. Rev. Lett. 93 (2004) 013401. 

8 P. Koskinen, H. Hakkinen, G. Seifert, S. Sanna, T. Frauen- 
heim, M. Moseler, Density-functional based tight-binding 
study of small gold clusters, New Journal of Physics 8 
(2006) 9. 

9 P. Koskinen, S. Malola, H. Hakkinen, Self-passivating edge 
reconstructions of graphene, Phys. Rev. Lett. 101 (2008) 
115502. 

E. Bitzek, P. Koskinen, F. Gahler, M. Moseler, P. Gumbsh, 



Structural relaxation made simple, Phys. Rev. Lett. 97 

(2006) 170201. 

11 S. Malola, H. Hakkinen, P. Koskinen, Raman spectra of 
single-walled carbon nanotubes with vacancies, Phys. Rev. 
B 77 (2008) 155412. 

12 C. Kohler, G. Seifert, T. Frauenheim, Density functional 
based calculations for Fe n (n < 32), Chem. Phys. 309 
(2005) 23. 

13 T. Frauenheim, G. Seifert, M. Elstner, T. Niehaus, 

C. Kohler, M. Amreutz, M. S. Z. Hajnal, A. D. Carlo, 
S. Suhai, Atomistic simulations of complex materials: 
ground-state and excited-state properties, J. Phys.: Con- 
dens. Matter 14 (2002) 3015-3047. 

14 G. Seifert, Tight-binding density functional theory: an ap- 
proximate kohn-sham DFT scheme, J. Chem. Phys. A 111 

(2007) 5609. 

15 N. Otte, M. Scholten, W. Thiel, Looking at self-consistent- 
charge density functional tight-binding from a semiempir- 
ical perspective, J. Phys. Chem. A 111 (2007) 5751. 

16 M. Elstner, SCC-DFTB: What is the proper degree of self- 
consistency?, J. Phys. Chem. A 111 (2007) 5614. 

17 G. Seifert, H. Eschrig, W. Bieger, Eine approximative 
variante des LCAO-Xa-verfahrens, Z. Phys. Chemie 267 
(1986) 529-539. 

18 W. M. C. Foulkes, R. Haydock, Tight-binding models and 
density-functional theory, Phys. Rev. B 39 (1989) 12520. 

19 O. F. Sankey, D. J. Niklewski, Ab initio multicenter tight- 
binding model for molecular-dynamics simulations and 
other applications in covalent systems, Phys. Rev. B 40 
(1989) 3979. 

20 T. Frauenheim, F. Weich, T. Kohler, S. Uhlmann, 

D. Porezag, G. Seifert, Density-functional-based construc- 
tion of transferable non-orthogonal tight-binding poten- 
tials for Si and SiH, Phys. Rev. B 52 (1995) 11492. 

21 G. Seifert, D. Porezag, T. Frauenheim, Calculations of 
molecules, clusters, and solids with a simplified LCAO- 



22 



DFT-LDA scheme, International Journal of Quantum 
Chemistry 58 (1996) 185-192. 

22 J. Widany, T. Frauenheim, T. Kohler, M. Sternberg, 
D. Porezag, G. Jungnickel, G. Seifert, Density-functional- 
based construction of transferable nonorthogonal tight- 
binding potentials for B, N, BN, BH and NH, Phys. Rev. 
B 53 (1996) 4443. 

23 F. Liu, Self-consistent tight-binding method, Phys. Rev. B 
52 (1995) 10677. 

24 T. N. Todorov, Time-dependent tight-binding, J. Phys.: 
Condens. Matter 13 (2001) 10125-10148. 

25 B. Torralva, T. A. Niehaus, M. Elstner, S. Suhai, 
T. Frauenheim, R. E. Allen, Response of Ceo and C n to 
ultrashort laser pulses, Phys. Rev. B 64 (2001) 153105. 

26 T. A. Niehaus, D. Heringer, B. Torralva, T. Frauenheim, 
Importance of electronic self-consistency in the TDDFT 
base treatment of nonadiabatic molecular dynamics, Eur. 
Phys. J. D 35 (2005) 467-477. 

27 J. R. Reimers, G. C. Solomon, A. Gagliardi, A. Bilic, N. S. 
Hush, T. Frauenheim, A. D. Carlo, A. Pecchia, The green's 
function density functional tight-binding (gDFTB) method 
for molecular electonic conduction, J. Phys. Chem. A 111 
(2007) 5692. 

28 T. A. Niehaus, S. Suhai, F. D. Sala, P. Lugli, M. Elst- 
ner, G. Seifert, T. Frauenheim, Tight-binding approach to 
time-dependent density-functional response theory, Phys. 
Rev. B 63 (2001) 085108. 

29 Hotbit wiki at https : //trac . cc .jyu.fi/projects/hotbit. 

30 http : //www . gnu. org/licenses/gpl . html. 

31 ASE wiki at https://wiki.fysik.dtu.dk/ase/. 

32 R. G. Parr, W. Yang, Density-Functional Theory of Atoms 
and Molecules, Oxford university press, 1994. 

33 N. Bernstein, M. J. Mehl, D. A. papaconstantopoulos, 
Nonorthogonal tight-binding model for germanium, Phys. 
Rev. B 66 (2002) 075212. 

34 R. S. Mulliken, J. Chem. Phys. 23 (1955) 1833. 

35 J. Junquera, O. Paz, D. Sanchez- Portal, E. Artacho, Nu- 
merical orbitals for linear scaling calculations, Phys. Rev. 
B 64 (2001) 235111. 

36 J. C. Slater, G. F. Roster, Simplified LCAO method for 
the periodic potential problem, Physical Review 94 (1954) 
1498. 

37 J. A. Pople, Molecular-orbital theory of diamagnetism: I. 
an approximate LCAO scheme, J. Chem. Phys. 37 (1962) 
53. 

38 T. B. Boykin, R. C. Bowen, G. Klimeck, Electromag- 
netic coupling and gauge invariance in the empirical tight- 
binding method, Phys. Rev. B 63 (2001) 245314. 

39 M. Graf, P. Vogl, Electromagnetic fields and dielectric re- 
sponse in empirical tight-binding theory, Phys. Rev. B 51 
(1995) 4940. 

40 I. N. Bronshtein, K. A. Semendyayev, G. Musiol, 
H. Muehlig, Handbook of Mathematics, Springer, 2004. 

41 http://www.dftb.org. 

42 J. J. Mortensen, L. B. Hansen, K. W. Jacobsen, Real- 
spage grid implementation of the projector augmented 
wave method, Phys. Rev. B 71 (2005) 035109. 

43 GPAW wiki at https : //wiki . f ysik . dtu . dk/gpaw. 

44 J. P. Perdew, K. Burke, M. Ernzerhof, Generalized gradi- 
ent approximation made simple, Phys. Rev. Lett. 77 (1996) 



3865. 

45 H. Rydberg, M. Dion, N. Jacobson, E. Scroder, 
P. Hyldgaard, S. I. Simak, D. C. Langreth, B. I. Lundqvist, 
Van der Waals density functional for layered structures, 
Phys. Rev. Lett. 91 (2003) 126402. 

46 T. A. Halgren, Representation of van der Waals (vdW) 
interactions in molecular mchanics force fields: Potential 
form, combination rules, and vdW parameters, J. Am. 
Chem. Soc. 114 (1992) 7827. 

47 M. Elstner, P. Hobza, T. Frauenheim, S. Suhai, E. Kaxiras, 
Hydrogen bonding and stacking interactions of nucleic acid 
base pairs: A density-functional-theory based treatment, 
J. Chem. Phys. 114 (2001) 5149. 

48 P. Koskinen, L. Sapienza, M. Manninen, Tight-binding 
model for spontaneous magnetism of quantum dot lattices, 
Physica Scripta 68 (2003) 74-78. 

49 N. W. Ashcroft, N. D. Mermin, Solid state physics, Saun- 
ders college publishing, 1976. 

50 H. J. Monkhorst, J. D. Pack, Special points for brillouin- 
zone integrations, Phys. Rev. B 13 (1976) 5188. 

51 D. Frenkel, B. Smit, Understanding molecular simulation. 
From Algorithms to Applications, Academic Press, 2002. 

52 S. Malola, H. Hakkinen, P. Koskinen, Effect of bending on 
raman-active vibration modes of carbonanotubes, Phys. 
Rev. B 78 (2008) 153409. 

53 C. P. Liu, J. W. Ding, Electronic structure of carbon nan- 
otori: the roles of curvature, hybridization, and disorder, 
J. Phys.:Condens. Matter 18 (2006) 4077. 

54 V. N. Popov, Curvature effects on the structural, electronic 
and optical properties of isolated single- walled carbon nan- 
otubes within a symmetry-adapted non-orthogonal tight- 
binding model, New J. Phys. 6 (2004) 17. 

55 P. W. Atkins, R. S. Friedman, Molecular Quantum Me- 
chanics, Oxford University Press, 2000. 

56 O. Kit, P. Koskinen, in preparation. 

57 I. Mayer, Simple Theorems, Proofs, and Derivations in 
Quantum Chemistry, Springer, 2003. 

58 B. Yoon, P. Koskinen, B. Huber, O. Kostko, B. von Is- 
sendorff, H. Hakkinen, M. Moseler, U. Landman, Size- 
dependent structural evolution and chemical reactivity of 
gold clusters, ChemPhysChem 8 (2006) 157-161. 

59 F. Yin, J. Akola, P. Koskinen, M. Manninen, R. E. Palmer, 
Bright beaches of nanoscale potassium islands on graphite 
in STM imaging, Phys. Rev. Lett. 102 (2009) 106102. 

60 A. J. Bridgeman, G. Cavigliasso, L. R. Ireland, J. Rothery, 
The mayer bond order as a tool in inorganic chemistry, J. 
Chem. Soc, Dalton Trans. 2001 (2001) 2095-2108. 

61 N. Bornsen, B. Meyer, O. Grotheer, M. Fahnle, E cov - a 
new tool for the analysis of electronic structure data in 
a chemical language, J. Phys.:Condens. Matter 11 (1999) 
L287. 

62 J. P. Perdew, Y. Wang, Accurate and simple analytic rep- 
resentation of the electron-gas correlation energy, Phys. 
Rev. B 45 (1992) 13244. 

63 D. D. Koelling, B. N. Harmon, A technique for relativistic 
spin-polarized calculations, J. Phys. C 10 (1977) 3107. 

64 R. M. Martin, Electronic structure Basic Theory and Prac- 
tical Methods, Cambridge University Press, 2004. 



